{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.012399,
     "end_time": "2021-04-05T17:55:52.635137",
     "exception": false,
     "start_time": "2021-04-05T17:55:52.622738",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Dosearch for Causal Identification in DAGs.\n",
    "\n",
    "\n",
    "This a simple notebook for teaching that illustrates capabilites of the \"dosearch\" package, which is a great tool. \n",
    "\n",
    "NB. In my experience, the commands are sensitive to syntax ( e.g. spacing when -> are used), so be careful when changing to other examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install rpy2\n",
    "%load_ext rpy2.ipython\n",
    "# import subprocess\n",
    "# subprocess.run('conda install -c conda-forge r-base', shell=True)\n",
    "# !pip install rpy2\n",
    "import numpy as np\n",
    "import statsmodels.formula.api as smf\n",
    "import rpy2\n",
    "from rpy2.robjects import r, pandas2ri\n",
    "pandas2ri.activate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "papermill": {
     "duration": 61.517411,
     "end_time": "2021-04-05T17:56:54.163802",
     "exception": false,
     "start_time": "2021-04-05T17:55:52.646391",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "install.packages(\"dosearch\")\n",
    "library(\"dosearch\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.012155,
     "end_time": "2021-04-05T17:56:54.189437",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.177282",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We start with the simplest graph, with the simplest example\n",
    "where $D$ is policy, $Y$ is outcomes, $X$ is a confounder:\n",
    "$$\n",
    "D\\to Y, \\quad X \\to (D,Y)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.013694,
     "end_time": "2021-04-05T17:56:54.215310",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.201616",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now suppose we want conditional average policy effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.202016,
     "end_time": "2021-04-05T17:56:54.429124",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.227108",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "graph <- \"x -> y\n",
    "  x -> d \n",
    "  d -> y\"\n",
    "dosearch(\"p(y,d,x)\", \"p(y | do(d),x)\", graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.012715,
     "end_time": "2021-04-05T17:56:54.454351",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.441636",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "This recovers the correct identification formula for law of the counterfactual $Y(d)$ induced by $do(D=d)$:\n",
    "$$\n",
    "p_{Y(d)|X}(y|x) := p(y|do(d),x) = p(y|d,x).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.043133,
     "end_time": "2021-04-05T17:56:54.510274",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.467141",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "dosearch(\"p(y,d,x)\", \"p(y | do(d))\", graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.01273,
     "end_time": "2021-04-05T17:56:54.535701",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.522971",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "This recovers the correct identification formula:\n",
    "$$\n",
    "p_{Y(d)}(y) := p(y: do(d)) = \\sum_{x}\\left(p(x)p(y|d,x)\\right) \n",
    "$$\n",
    "We integreate out $x$ in the previous formula.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.013627,
     "end_time": "2021-04-05T17:56:54.562271",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.548644",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Suppose we don't observe the confounder. The effect is generally not identified.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.046266,
     "end_time": "2021-04-05T17:56:54.621562",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.575296",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "data <- \"p(y,d)\"\n",
    "query <- \"p(y | do(d))\"\n",
    "graph <- \"x -> y\n",
    "  x -> d \n",
    "  d -> y\"\n",
    "dosearch(data, query, graph)"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAEwCAYAAAAQDc26AAAgAElEQVR4Ae2dB/QURbrFQST8CYoiOQeVRYFFQFBXQRdFySoKuqiAWaIIghnf+gwE9QlKkKCLuiuCKKhPFFlQESUJC2tgJUpGJUhGrXduv63e7pnpme6eDtXdt86Z0zMdKvyqq+d21VdfFRIMJEACJEACJEACJEACgRIoFGhqTIwESIAESIAESIAESEBQgPEmIAESIAESIAESIIGACVCABQycyZEACZAACZAACZAABRjvARIgARIgARIgARIImAAFWMDAmRwJkAAJkAAJkAAJUIDxHiABEiABEiABEiCBgAlQgAUMnMmRAAmQAAmQAAmQAAUY7wESIAESIAESIAESCJgABVjAwJkcCZAACZAACZAACVCA8R4gARIgARIgARIggYAJUIAFDJzJkQAJkAAJkAAJkAAFGO8BEiABEiABEiABEgiYAAVYwMCZHAmQAAmQAAmQAAlQgPEeIAESIAESIAESIIGACSgtwE4++WRRUFAgihcvLooWLSqKFCkiTjjhBHHiiSdq+0qWLCnKlCkjypYtK8qVKycqVqwoqlSpIqpXry5q1aolatSoof2uUKGCOOWUU0Tp0qVFiRIltM+yZcsCRs3kSCCaBH788Uet3aEdor2hrVWqVElrX2hjlStX1vaddNJJWttCO/3mm2/SCvvbb7+JOnXqaOfg3PLly4tq1aqJmjVr6nGgjaKt9+nTJ+167iABEiCBOBFQVoAdO3ZMFCpUyJcPBNzKlSvjVI8sCwn4RmD//v2iWbNm2stM4cKFLdskXo7q1q0rzj33XLFp06a0/Pz666+iRYsWol69etrLVKb2jfjr168vRo0alXY9d5AACZBAnAgoK8D27t2rP+jRo9WuXTvtrXj06NHi5ZdfFjNnzhRz5swR7777rnj//fczfmbNmqXHYXzYP/vss3GqQ5aFBAIjsHv3bjFo0KC0doVesZ9//tl2Po4cOSLOOOMMUzw33HCD2LFjh+04eCIJkIA1gaNHj4onnnhCDB8+XMyePVvgBcjrsGvXLq+jTFR8ygqwzZs36w/n5cuXu6oUPNCNwgvfe/To4SouXkQCJPAfAl27djW1LfRcZer1+s8V5m+ffvqpqRese/fuAkOUDCRAAvkT+Ne//iWaNm1qaqMNGzYUGzduzDty9Ij/+c9/FqeffrpmepB3hAmOQFkBBhstCCbYdbl5MD/33HOmmw9xNWnSRBw6dCjB1c2ik4A3BCC2YE9pfMHp16+frci3b98uqlatql/bqlUrgbd1BhIggfwJvPrqq5qtprFtyu8YTUIvttuAUSfEYYzPbVy8TghlBRgqGpV84403Oq6njz/+WDPUlzcJthgi8UL9O84MLyCBmBIYPHiw/iBGGytWrFjONoahx5YtW+rXNWjQQOzZsyemhFgsEgieACayGP/7Ur9fd911rjL12GOPpcWLXjYG9wSUFWATJkzQKhtbJ2HLli1ar5nxpsOsrHnz5jmJhueSAAnkIIDZkZiBbGxrvXv3znoVXqjk+Zg9yZeirLh4kAQcE3jyySfF3Llzxddffy1g74zZy7LNye2iRYscxTt06FBTHJhwc8EFF4hJkyY5iocnmwkoK8A++eQTcf/994v169ebc5zlF4YxMMtK3mRyO3LkyCxX8RAJkIBbAk899ZSpvVm5oED8mNko2yTe0t3adrrNK68jgSQS+Pvf/54mwjp16mQbBQSdbLfYtm3bVqxZs8b29TzRmoCyAsw6y9ZHbrvtNtONgpulW7du1hfwCAlEjMDBgweFSj7sYFMJX17GB/S1116bRvW9997TfInhPLiBwW8GEkgSAUws27BhQyhFfvDBB01tFJNmYKifK8DDAHq7ZPt+9NFHc13C4w4IxEaATZw4Ub9J5M2CWR/4w2IggbgQ6Ny5s3afjxs3TpkiTZ482dT28HD/8ssv9fzBKSucKst2ibbKQAJJIvDTTz9pzosxcQVDg0EHvCjBOblsg9g+8sgjWbPx3XffaQ7M5TUvvPBC1vN50DmBWAiwzz//PK2LFZ7v161b55wIryABhQnAiSkeiPAe/8svvyiRU+TjrLPOMj3c27dvr+UNBvZGf18PPPCAEnlmJkggSALobZJC5pZbbgkyaT2tVM8AeJZYBfgMO++88/Q8wwaMwXsCkRdgO3fuNE1px02OLlN0nTKQQNwI3HffffpDUaVeMDh6lH8wcrtw4UJx2WWX6fvhl4+BBJJKAG6Q0DZgFP/tt98GjmHfvn2iVKlSentEXpYuXZoxH3B4Ltsxet3duILKGDF3mghEWoAdP35cXHTRRfqNIm+Yxx9/3FRI/iCBuBDYtm2btqYp7nWsp5hpzcWwynrhhRemtUXZJi+55BKB5cUYSCCpBKZPn663DyzXBZcsQYc77rhDzwPa5sMPP5yWhbVr12prMOM4fH5h+JTBHwKRFmADBgww3Uy4Ya6++mp/SDFWElCEwPPPP6/f91h7Ea5XVAifffaZni8pvLA9++yzBZYWYyCBJBNAL1LHjh31NoLVJNCJEGRYvHixnj7aZqofL+RRvkhhRjN8ajL4RyCyAmzatGmmGwk3E5w6OlmPzj+sjJkE/COAh2TPnj31+x8iLIwhjUwlvPLKK/V8SRFGW8xMpLgviQTQm4TF5mXbwPBekBPF8OwwzlrGhBn0qsswZswYPW8PPfSQ3M2tTwQiKcAww6qgoEC/UXAzY5YVuk4ZSCAJBDCcB3cP8kGOSSewwwo7wPZS5klu4QwynwAj/ylTpggMn8D5I9agg+kBfmNNSQYSiBKB77//XruHZfv4/e9/b8slhFdl7N+/v6mNTp06VYsaazxKx8pnnnkmlwfzCniWeCInwOB9O3U6LVQ8li7yIrz55psCS6zA2Rwe9K1bt9Ye9PD469eYPQwh4SyPHzJwcg989NFH4qqrrtIfpmgHd955Z2i9wHiAw/WL/GOR21NPPdX1ckPoPTPOxpJxGreDBg3youmb4sC0fSd1wXPZdp3cA7NmzTLNHIZjYrhnQQ+V3+HDDz80tdEePXpoSY4YMULfz5Vj/K6F/48/UgIMU2MvvfRS/SaRD2EvnMPt2rVLwDuwjDPTtkaNGgIuL7wOjRs3zppuprxwXyEyK5SZAd5e8ZYdZEAvFVxPWN2X99xzj+Ps4IUHEw1knCVLltTaP2xU4MxV7scWq2Z4GeADyRg/v2e+18jFWy7du3f3fbIKVowpU6aMfn9XrVpV61yoUqWKts/tWpFetr+kxBUpATZs2DD9ppENH6LJi7cG3IQyTmxhgGj8Lb+ji9aOB2EnNxAFmLcPMVlXSd7COWqQYeDAgXp7geuJ4cOH679RD1io28myYsg7hhdxbcWKFcXMmTNNQyIwQ6hTp46eBtqvl4ECjG0yrOfH1q1bvbyVM8aVaquJUR+Ut2jRoqF568+Y0ZjvjIwAmzFjhv6wlQ0Db/rwbZJv2L17tx43FgtGV/Lhw4e1uN9++23Tgx5p51pw2Gl+MCUZHpL5IQM390DqywIM9DGEFlQYP3683n5gHvDDDz9ow6AVKlTQ96PdZFqiKFse5Rp0X3zxRcbTsB/DrvJ54KVLDgx9uqkLXsM2bPcegD8w4/2L+zgob/PwISjbjXF78803Z2xr3OkPgUgIsH/+85+67yN5s6AL1aslHWCDhXgbNWqUsfsXU+iNPWQYFmEggbAJYPaSdO6I+xeCJ2hDfNiTyOFA/PEYF9iG8b1sr3ILVxV2Q4cOHdKmyadea5xRtmrVqtTD/E0CShKYP3++bvCOttGiRYtAZzIbPfPLtokXOc5YDvZ2UV6AQfwYlzLBzYK3BhgxehXeeOMN7Y9iwYIFllFiFpa8UbHFWz4DCYRFAL6/5LJEuB8xO9A4nTyIfOEFSM6aQh5eeuklU7Kw4YLdpLHdwKDebkAbxx9VtgAfY4j/tNNOC7TXL1ueeIwEshHATGG8rMh2gVmJYTgpNrqjQF7atGmTLds85gMBpQUYbLuMjuvkDYuV3b0M6I7FjMdstmTohZPpY6uK80svOTCuaBCArzv01sr7EQ4d/Zqha0UELyDwPybzcNddd2U8NfXFBefDI7gXAcuQyd63fF1deJEfxkECuQisWLFCXw4IS+bB71ZYoVWrVnr7RbsMavgzrPKqmK7SAgyzG+UDXm7btWsnMBvSTcAiqDAWdhPQNSvzgIe+KgshuykLr4k2AdhpyHsRjhyDvhcxi8q4BBh6taze4JE34zAh8g3jecSRb+jWrZvG4fLLL/ckvnzzw+tJIBsBOFzFi75su2PHjs12uu/HUnunVXHm7HvBFUpAWQH2zjvvpBkoYshlz549rvChxwpDlxiycBOQH9lw4AqDgQTCIAB7RWm4C2Fz4MCBwLNh9MKPGYq5Zm3JIX7ZfrDFYr9uA2ZTdunSRWuPMOy3En9u4+d1JOAHgccee0z/D7n11lv9SMJ2nPDZJ58jaI+Y/Rj0ski2MxvjE5UUYPBoD8/2xgc2VnFfs2aN66qQSyw4nYklE8QakzI/b731ltzNLQkESgCzdOV9mM1m0a9MyZmJyAN6gu2sFYehfaw5J/ONLTz3w6my3YDZyOXLlxflypUzxVOzZk2BPzOugmGXJM8LgwB6gitXrqzdu3hpgQAKM2AGsbE9/u53vwszO4lNWzkBhjf6s846y3Rz4EbBW3Q+QS4w6sZpK6a3y6n+8J/CQAJhEYAIQXtwYszuVV6xSoTxrdmJ3dXcuXPT2rQTc4B7771XE3wwFL7++utF8+bNBWxo5J8I3uDhPoaBBFQksGzZMv1edfMf5HWZUm0z0aPMEDwB5QTYNddco9+o8uE6dOjQvMhgyq3843Aq5PDm0rJlSy1P1atXFzD8ZSCBsAg8/fTTAm+rX331VaBZwNAnPNHLNummJ/n888/Xr0c8EE127U5gg5na9tD7Jr13Iz70qm3fvj1QLkyMBOwQgF8+vLxjJCXs3i/kFytTyLaMLZycMwRPQCkBZlyLSt4cGHrEg3bRokW2Ppi2Dl9Ir7zyisCY+w033GAyfMSbiJMgPXpjSDSfIVAnafJcElCFAOxC0JawVp1sk9g+88wzjrOYyS8Y7DrzWTYJriqM+Xrttdcc54sXkIBfBDBhDI6+VQtYrcLYblJdyKiW37jmRxkBtnjxYn2Yz3hjeP3dia8kTJdHzxk8FufyRxTXG4TlSiaBTz75RGDGMXqVrNogFtnGVPa+fftaQsLMYwwbnnPOOWkiTsYLWzLYiGEYxI5NWWpiRifJ9913X+ph/iaB0AhcccUV2lD5qFGjQstDpoRhhybbH7b4/2UInoAyAgwK3HhD+PEdQsquCwusQwdnefhzoNF98DcmUwyXgJP2iKF5q5BNwGVq427amtGf0UMPPWSVFe4ngcAJ1K5dW/9fg09LFSaLbNiwQc+TbINOJsQEDjHGCSojwOBccuPGjZqDU9h64IbAOo8wyse6jJhqDnssCCjMqsIHv7EfPoVwDsbZcT7G2H/66Set63fHjh1anJs2bbI96wpDjXi7R+8Xhl8YSCBpBGR7RI8xnK6iTcHZK9oc2hpWqNi1a5c2fLh582ZLPGjTGGJEm0abRPtEe5XxoI1jiAbpGOOBw8p58+ZZxms8AIN8+UeCRbsZSEAVAuhZMr6E4IX+qquusn1v+1EOufC2bDPYrl692o+kGGcOAsoIsBz5DOwwjH2lYS89AweGnQmRgIkAfI3BUSSEWraAFzGszYo/ESxH5NZPYLY0eIwE8iGAlxA4TDYKHnzHjOa2bdtqBvAwd/nHP/7hi18/uJzApBnYfRkdwRrzA2GIyWYwA1BtuDQf9qpfSwFmqCE4lISXbtyY8HfEQAIkEA4B/CGhHeYyDjYOlU6aNCmczDJVErBBAAvRYxakdGlkFEDye4MGDTxf09XYRmQ62bZw88IQDAEKsH9zxjALbn7cmPfff78lfQyBYqiTb9qWiHiABPImgKFJ/FFhFjRmQGcKcC9TqVIlrc327t0761quma7nPhIIgwCG27H+cOPGjXX3SEZBtGTJEk+zhV5i9MIhXQz3o23hPwwznGHCA7MAmAfAXAemOvgvZAiGAAWYEJp9S7NmzbQHORpCkyZNNDGGxYYxwwpDG5iGD79FxoYCQ30GEiABfwiceeaZWnuDw9Xbb79dm4mMJcXgSmbkyJFam8QxOGnFnwwDCahOAPaTzz33nG7mYvw/wXc6+la9Br3NX+IFGBqEcRZVaoPI9tsrAYYJB/n4QvL2lmBsJKAGAfRqWbU/TJBp3bq18Lq3QI2SMxdxI5BLeOE+L1OmTJqz4bhxYHnMBBIvwOCJW3rJt3rYW+2HryQvAjz9w0VGnz59KMS8AMo4YkMAMyOnTp0qHnzwQQGfYvDg/fzzzwssyM1AAqoTyCa80HuLJcXk/wtdqKhem97nL/ECDEgx5o0P7Low/R7j4xgbh70XPviOfTiG6fc4F1PwvQiIx+hlnELMC6qMgwRIgATCI5BLeHXt2lWb9di+fXtNgBUrVoy9X+FVV2gpU4CFhv7/E8YixXIavXwTwpZCLOSKYfIkQAIk4JCAXeGFaPFCD+GF5z3cVDAkjwAFmAJ1DhuwBx54gEJMgbpgFkiABEjAKQEpvIzLYskXagw1yh4vY7wffPCBPvw4ceJE4yF+TwgBCjCFKppCTKHKYFZIgARIIAcBN8JLRjl69GhdgK1atUru5jZBBCjAFKxsO0IM0/EZSIAESIAEgieQj/CSuR00aJAuwA4ePCh3c5sgAhRgCld2NiGGhcL79u2rrXOpcBGYNRIgARKIDQEvhJeEceedd2oCDHZgDMkkQAEWgXqnEItAJTGLJEACsSXgpfCSkO6++25NgGEdxlxrnspruI0XAQqwCNUn3F9gmaRMsybZIxahimRWSYAEIkHAD+ElCz5+/HhNgP3ud7/T3B3J/dwmhwAFWATrWgoxeE6WM23klkIsghXKLJMACShFAMJrzJgx2lJ08tkqt1azGp0WAD4mP/roI22dRqfX8vx4EKAAi3A9UohFuPKYdRIgAeUIBCG8lCs0MxQaAQqw0NB7lzCFmHcsGRMJkEDyCFB4Ja/OVSgxBZgKteBRHijEPALJaEiABBJB4OjRo1mHGq+++mqRBB9dy5YtExdffLFYuHBhIupdlUJSgKlSEx7mg0LMQ5iMigRIIHYEILzGjRsnatSokWZHCxuvpAgvVCzWQJb2bdhu3bo1dvWtaoEowFStGQ/yRSHmAURGQQIkEBsC2YRX4cKFEyW8ZKX27t3bJMA6deokD3HrMwEKMJ8BqxA9hZgKtcA8kAAJhEWAwisz+RdffNEkvmRP2JNPPpn5Au71lAAFmKc41Y4slxDr168fPeurXYXMHQmQgAMCuYRX586dE2HjlQnZjBkzBJzAStFl3GIYdtKkSZku4z4PCVCAeQgzKlFRiEWlpphPEiABNwTsCK/ly5e7iToW18yZM0cULVpUE19FihQxibDixYtrvyHCpk2bFovyqloICjBVayaAfEGI3XfffcLKoSt7xAKoBCZBAiTgGQEKr+wojx8/Lv7rv/7LJL6kR37ZAzZhwgSB9SnxG+JsyJAh4siRI9kj5lFXBCjAXGGL10UUYvGqT5aGBJJGAMILQiLTrEYY12OoMck9XvC6/7e//U3Ur19f7+0CF4gtrDUsxRe2a9euFdOnT9fEl9xfu3ZtMWXKFAEBx+AdAQow71hGPiYKschXIQtAAokiQOGVu7ohmrp162YSWZUrVxbvvvuudnEmAYYD8AkG4SVFGLaXX365OHToUO5EeYYtAhRgtjAl6yQ7Qoy+YpJ1T7C0JKASAQov+7WB57k0ti9btqw2BLlv3z49AisBhhMOHjwoRo0aJSpUqKALsXXr1unX8kt+BCjA8uMX66spxGJdvSwcCUSOQC7hBR9WSR5qtKrQlStXCric2L9/f9op2QSYPBlCbOrUqWLx4sVyF7ceEKAA8wBi3KOgEIt7DbN8JKA2AQov/+rHjgDzL/Vkx0wBluz6d1R6CjFHuHgyCZBAngQovPIEaONyCjAbkHw6hQLMJ7BxjpZCLM61y7KRQPgEpPCqWbOmbnskjcExe49Djd7VEQWYdyydxkQB5pQYz9cJ5BJi/fv358KuOi1+IQESyEWAwisXIe+PU4B5z9RujBRgdknxPEsCu3fvzurQlULMEh0PkAAJCCEovMK7DSjAwmNPARYe+9ilTCEWuyplgUjAVwIUXr7itRU5BZgtTL6cRAHmC9ZkRwohNmzYMMsljtgjluz7g6UnAQovde4BCrDw6oICLDz2sU+ZQiz2VcwCkoAjAhRejnAFcjIFWCCYMyZCAZYRC3d6SUAKsdKlS6fNaCpRooRgj5iXtBkXCahHgMJLvTqROaIAkySC31KABc88sSlSiCW26lnwhBKA8MKCz9ncSSxbtiyhdNQoNgVYePVAARYe+8SmTCGW2KpnwRNCgMIrOhVNARZeXVGAhcc+8SlTiCX+FiCAmBGg8IpehVKAhVdnFGDhsWfK/yaQTYgVFBTQRox3CgkoTiCX8OrYsaPgUKOalUgBFl69UICFx54ppxDIJcQGDBhAz/opzPiTBMIkQOEVJn1v0qYA84ajm1gowNxQ4zW+ErAjxLZt2+ZrHhg5CZCANQEKL2s2UTtCARZejVGAhceeKecgQCGWAxAPk0DABCi8AgYeQHIUYAFAtkiCAswCDHerQwBCbOjQoSKTHzHYiGFokj1i6tQXcxI/AhRe8atTWSIKMEki+C0FWPDMmaJLAhRiLsHxMhJwSYDCyyW4CF1GARZeZVGAhceeKbskQCHmEhwvIwGbBCi8bIKKwWkUYOFVIgVYeOyZcp4EKMTyBMjLSSCFAIVXCpAE/KQAC6+SKcDCY8+UPSJAIeYRSEaTWAIUXomtekEBFl7dU4CFx54pe0yAQsxjoIwu9gSk8KpVq5YoVKiQ6VO4cGFBB6qxvwUowEKsYgqwEOEzaX8I7Nq1i7Mm/UHLWGNCAMJr4sSJgsIrJhWaRzHYA5YHvDwvpQDLEyAvV5dALiE2cOBAuq9Qt/qYMx8I2BFeS5cu9SFlRqkqAQqw8GqGAiw89kw5IAIUYgGBZjLKEqDwUrZqQs8YBVh4VUABFh57phwwAQqxgIEzudAJHDt2LOtQY4cOHQR7vEKvplAzQAEWHn4KsPDYM+WQCECI3XvvvZae9Tk0GVLFMFnPCFB4eYYy9hFRgIVXxRRg4bFnyiEToBALuQKYvOcEKLw8Rxr7CCnAwqtiCrDw2DNlRQhQiClSEcyGawIUXq7RJf5CCrDwbgEKsPDYM2XFCFCIKVYhzE5OAhReORHxhBwEKMByAPLxMAWYj3AZdTQJ2BFi27dvj2bhmOtYEKDwikU1KlEICrDwqoECLDz2TFlxAhRiildQArNH4ZXASve5yBRgPgPOEj0FWBY4PEQCICCFWKlSpUxLtWDploKCAoFZk+wR473iJwEKLz/pJjtuCrDw6p8CLDz2TDliBCjEIlZhMcguhVcMKlHxIlCAhVdBFGDhsWfKESWQS4jdfffd7BGLaN2qkm0KL1VqIv75oAALr44pwMJjz5QjToBCLOIVqGD2IbxefPFFy0Wy6blewUqLeJYowMKrQAqw8Ngz5ZgQgBAbMmSIsLIRY49YTCrax2JQePkIl1FnJUABlhWPrwcpwHzFy8iTRIBCLEm17U1ZKby84chY3BOgAHPPLt8rKcDyJcjrSSCFAIVYChD+TCOQS3i1b99eLFmyJO067iABrwlQgHlN1H58FGD2WfFMEnBEIJsQK1mypPBqaPL48ePiwIEDjvLGk8MhIIVX7dq101yaFC5cWFB4hVMvSU6VAiy82qcAC489U04IgZ07d1raiOUrxPbu3SuaN28u8IeOdBjUJEDhpWa9MFdCUICFdxdQgIXHniknjIAfQmz+/Pl6T8rs2bMTRlT94lJ4qV9HSc/hnj17RI0aNfTP+vXrk44ksPJTgAWGmgmRwP8T8FKIQXTBIz8+EGMMahCg8FKjHpgLElCZAAWYyrXDvMWagBdC7L333tMF2Jw5c2LNKwqFyya8IJJp4xWFWmQeSSAYAhRgwXBmKiRgSSCXEBs0aJClZ/3ly5frAmzs2LGWafCAvwQovPzly9hJII4EKMDiWKssUyQJuBFiP//8szjhhBM0Eda7d+9IljvKmabwinLtMe8kEC4BCrBw+TN1Ekgj4FSINWnSRBNg1atXF7/99ltafNzhPQEKL++ZMkYSSBoBCrCk1TjLGxkCEGKDBw/OuMQR3FdgaHLHjh3ioYce0ochFy5cGJnyRTGjFF5RrDXmmQTUJEABpma9MFckoBPIJcR69eol4MQTRt5du3bVr+MX7whAeE2aNEnztyZnnRq3NK73jjVjIoGkEKAAS0pNs5yRJ5BNiBUpUkQTYLAHW7VqVeTLqkoBKLxUqQnmgwTiR4ACLH51yhLFnEA2IYZemcaNG9MWLM97gMIrT4C8nARIICcBCrCciHgCCahJYMWKFaJjx46iaNGiug2YHBYbN26cmpmOSK6uv/76NKZg265dOy6SHZE6ZDZJQHUCFGCq1xDzRwIGAr/++quA9/vWrVtnFAhSgD377LOGq/jVKQHjEk9SeH3xxRdOo+H5JEACJGBJgALMEg0PkIBaBD744APRqFGjNOFVqlQpce655woY448YMULMmjWLQ5AeVF2rVq20Hi8KLw9gMgoSIIE0AhRgaUi4gwTUIvDLL7+I/v37m4RXQUGB6NGjh5g7d67AcT/Djz/+qDl7LV68uChTpowoV66cqFSpkr54b+XKlbV9J510kihRooTAhIBvvvkmLUvwUVanTh3tHJxbvnx5Ua1aNVGzZk0h4yhdurQ2pNqnT5+064PecejQoaCTZHokQAIJIkABlqDKZlGjSeDhhx/WxReEF/x/wRA/qLB//37RrFkzUaVKFRfzczkAACAASURBVN3dhRzqNG4xA7Nu3bpab9ymTZvSsofh0xYtWoh69eppIs14rfwOdxr169cXo0aNSrueO0iABEggTgQowOJUmyxLLAncc889mgCDp/vVq1eHWsbdu3drAlAKJrlFrxiWRbIbjhw5Is444wxdWCKeG264QXMsazcOnkcCJEACUSZAARbl2mPeE0Fg3759YsaMGQJDgaoEOHyV4gtb9Fxl6vWyyu+nn35q6gXr3r077dasYHE/Cdgg0LRpU214/5RTTtF6q/HCdtppp2kracAsYMGCBTZisX/K008/rbVhmA3AnADpwZSgbNmyAuYKWKGDITsBCrDsfHiUBEggAwGILdh7GUVYv379MpyZvmv79u2iatWq+rUwdj969Gj6idxDAiRgm0Dnzp1FrVq1xIknnqi3LWP7hCmDVwE92DBJMMYvv5966qkC69NyJnZu2hRguRnxDBIggQwEsE6lfOhiW6xYMbFx48YMZ/5nFx7cLVu21K9r0KCB2LNnz39O4DcSIIG8COBl5vXXX09bQ/aCCy7IK17jxbDRNLZ9+R0zsGHryWCPAAWYPU48iwRIIIUAhkQx3CAfvtj27t075SzzzxtvvFE/H8MVuQSb+Wr+IgESsEvgoosu0tsa2iYcNh84cMDu5Zbn/fDDDwLDnMZ2j+9YgYPBGQEKMGe8eDYJkICBwFNPPWV6EFu5oMAlxrdm2I0sX77cEBO/kgAJeEUArmkwFJgqkv73f/837yT69u2bFi/SufTSS/OOO2kRUIAlrcZZXhLwkAB8ZcGXl/FBf+2116al8N5772m+xHAebFTwm4EESMAfAp999pmpTcr2CbOBfMKaNWssbcyuu+66fKJO5LUUYImsdhaaBLwjMHnyZNPDHjMiv/zySz0BOGU9+eST9XMmTpyoH+MXEiAB7wkYfQdK8YXtOeeck1dibdq00dsx7DeNcdudhJNXBmJ2MQVYzCqUxSGBoAlguOOss84yPYzbt2+vZQMG9kZ/Xw888EDQ2WN6JJA4As2bN9faI1auwEcKJThLduvOBgb2Mp4LL7xQSP+Ect+jjz6aOM75FpgCLF+CvJ4ESEBbIFw+iOV24cKF4rLLLtMf2nC0ykACJOAvAaySgV5otMMuXbqIK6+8Um+D2Ddz5kzHGTh8+LCoXbu2Fg/sPFetWqVNuJFtHduxY8c6jjfpF1CAJf0OYPlJwCMCeCs2PpCN3y+55BJx7Ngxj1JiNCRAAlYEpk2bprdDiKIxY8bov9Em3ayz+thjj+lxwAgfAeLO2Mb/9re/WWWJ+y0IUIBZgOFuEiABZwSsDH/PPvtssXfvXmeR8WwSIAFXBK6//npdGMH+EobzRqGEtVadhC1btug+xTCzUg5hpr5wzZs3z0m0PFcIQQHG24AESMAzAqnDHXjwr1u3zrP4GREJkIA1AThBxfJDaHdYGkiGihUrmkTYtm3b5KGcW6OgQ2+aDKl2nytXrpSHuLVJgALMJiieRgIkkJvA+++/b3rQ44+AS5Lk5sYzSMALAp9//rne/nr16qVHibVWjb1gr7zyin4s2xes2Sqvg+DChBsZ4EhZHsP2+++/l4e4tUmAAswmKJ5GAiSQncD+/ftFw4YNTQ9lPJgxbJHPckN4sI8YMULAiz6m0WNIE4uBY7FfvnVnrxMeTRaBRx55RG9/r732ml54uH4xiiWjONNPSvmC3jS0N3ndhx9+aDoDC27LY9jCUJ/BGQEKMGe8eDYJkEAGAngzhusJ4wPZ+B1T1t0EGPZmWvZExo3ZXrfffjvXn3MDl9fEjkCLFi20Noh2gdmQMnz33XemtolFu3MFo2jDQt/GcPDgQVN8pUqVMh7md5sEKMBsguJpJEAC1gQGDhyoP5DhemL48OH6b4glLNS9fv166wgyHHnxxRdNcSAeeNGX4su4lTOzMkTDXSSQCAK7d+/WV5to1KhRWplr1KhhajsQZVYBPdbly5fXzkfbTT138+bNprhq1qxpFRX3ZyFAAZYFDg+RAAnkJjB+/Hj9YYw3ayzW+/PPP4sKFSro+yGWMi1RlC12TJfHdTAgnjp1qtiwYYPW04U/gzvvvFP/s8E58E3kVOBlS5vHSCBqBF599VW9vWXqce7Zs6d+HG0m24oUxheqe++9Nw0Fhv4Rh/zk62E/LYGE7KAAS0hFs5gk4AcB2IXIXqkSJUqYFtiG8b18QMstXFXYDXJI8/XXX894yciRI03xP/300xnP404SSAKBHj166O0h06Lbf/nLX/TjaI8wzM8UvvrqK71NV6pUScC2MzXMnz/fFBd6vRmcE6AAc86MV5AACQghvv76a1G2bFn9QfzSSy+ZuBw5ckSkDnucd955pnOy/cCsq9atW1ueArsz6Z0bfygchrRExQMxJwCDeeOQIWy0UgMms8gXIWzRs5wpGFevmDJlSqZTxIwZM0xxwVUFg3MCFGDOmfEKEkg8AQwz1q1bV38I33XXXRmZ4AFufOjj+/Tp0zOem7oT4mrChAmpu02/r7nmGj3+W265xXSMP0ggKQSWLFmit4NsLy3GdVnRFlevXm1C9Pbbb+vxNG3aVPz222+m4/KH0UAf8fTv318e4tYBAQowB7B4KgmQgBBHjx4VF110kf6gRq+W1TJD6KWC522jCKtTp44Whxcssb6kjJuLAXtBlHFEkQDufdkOsGyQVbjjjjv083D+//zP/+inosfa+FIFH2BW4cknnzTFw7ZnRSr7fgqw7Hx4lARIIIWA0ZgXwxhbt25NOcP884033jA9rPHgHz16tPkkl7+aNWumx71o0SKXsfAyEog2gZYtW+rt4IsvvrAsDHqfpVDDtlOnTvq5jz/+uH7suuuu0/dn+jJkyBD9XMTz/PPPZzqN+3IQoADLAYiHSYAE/kPA+OYL4/uPP/74PwctvmEYA8MZxgc/fHvJNeUsLsu5+x//+IceZ5MmTXKezxNIII4EYA5wwgknaG0BNplGb/Wp5d21a5eAjzDZFuX5eIkqXbq0tr9kyZICbiayhZtvvlmPA3FZTZTJFgePcS1I3gMkQAI2Cbz55pumh7eTJYbmzp1remDjoY2p7vkEOesLQnDZsmX5RMVrSSCyBP7617/qbQtrseYK8BEmBRi26DGTbQm/4cMvV0hd8/Wjjz7KdQmPZyDAHrAMULiLBEjATGDp0qUCb8bywe3UpxdiO//88/XrEU/RokXFt99+a07I5q85c+bocT3xxBM2r+JpJBA/AliiS7bLF154IWcBjT6+cF3Hjh31FyvMWj506FDOOFq1aqWniTi4JFhOZBlPoADLiIU7SYAEQOD48eMCC/fK4Qn5oH/mmWccA8rkF6xevXqOF/HFEiuwPUNeOPPRcTXwghgRwPC+bAtoD2vXrs1ZutmzZ5vEk2zT2GLpLzshdc3XLVu22LmM56QQoABLAcKfJEACQnzyySeiXbt2WddhxCLbeBPO5n8LAqlNmzbaor6pIk4++DGECBuxLl265LQpg3+jc889V/sDQf4gEBlIIKkEFi9erIsp9F7ZCXv37tVWjpDtT27/8Ic/2LlcO6dKlSp6urieC3HbRmc6kQLMhIM/SIAEQABOVeWDOde2evXqltCyLaSdKd633nrLMi44m8SiwLgObjDsDJVYRsYDJBADAliIXrYjO/ZfssjNmzfXr8P1MMy3a0cJI//ixYubrueLkCTrbEsB5owXzyaBRBDAWo4bN24U27Zt09Z2xHIk8BOEhy/edvEWjRlV8K6dbcYU4sA5GDb86aefxIEDBzQfYDKeffv2CSwijHSyxQPo0ocRXE/gOgYSSDIBOE3FGqhSgDVu3Fhbg9UOk6FDh+rX4Xq7nuwhtMaMGWO6FtePGDGCvdF2wKecQwGWAoQ/SYAE1CMwbNgw7aHfoEEDTRCql0PmiASCIYAXkVS3LlKEYTgfvVtXXXWVwMQZq2CclYzeLCx0bxXQM4aeZyy4XaZMmTTxJdOGicHFF18s/vSnP1lFxf0pBCjAUoDwJwmQgFoEMMsRD3l40EdPGQMJJJUAbCCNs5Gl+Mm0nTdvniUmxFOsWDGtXQ0aNMjyPBxIXfcxU1rGfVjAm8EeAQowe5x4FgmERgCGtnij/eCDD0LLQ1gJY1o9Hu5Vq1YV69evt8wGlkLCMCmNgS0R8UBMCGDIHo5T4YAVpgJYGgz3P4b34dwYLyno0cL+bAErWsA1jB2HyJs2bdLiRdoY/kc7gxmBTBfmBTt27NDMDWCewGCPAAWYPU48iwRCIyCX/oG9x/3332+57mJoGfQpYbi/kB6+MbwBB5JYTLhmzZra1Ht48S4oKNDPgVDDjEsGEiABEogCAQqwKNQS85hoArDXgNNS2c1fu3ZtAT9c6PGJa4CvItizyDLb3VKAxfWOYLlIIH4EKMDiV6csUQwJwC9XrVq1TIIEvUIdOnTQxNjq1asFnDLGJWBWll3RZTzvj3/8o2cIFixYkHOhcc8SY0QkQAKJI0ABlrgqZ4GjSgCGsyNHjhQVKlTIKE4gyMqXLy+uu+46sX379qgWU8s3bEzg5gK2JbA5QdnhBgPT4CE0scVv7MdxnAf7FPz2IsDHGJxNlihRQvTv359CzAuojIMESMBEgALMhIM/SEB9AhAHU6dOFS1atMgoxNAjBNHA4J5A6rJJFGLuWfJKEiCBzAQowDJz4V4SUJoAZjiNGzdOVKtWLaMIy+YDSOmCKZI59DSiR9E4vInvMPpnj5gilcRskEDECVCARbwCmf1kEYDwGj9+vMC6b6niQP6+9957kwXFp9JiSBMOYCnEfALMaEkg4QQowBJ+A7D40SCQTXhhHbdSpUppgqxcuXJcpsfjKs0lxAYMGEAbMY+ZMzoSSAIBCrAk1DLLGFkCuYRXp06dxHPPPaf3hj388MORLavqGbcjxOipX/VaZP5IQB0CFGDq1AVzQgI6ATvCa/ny5dr5N910kybA4LQ014LWegL84poAhZhrdLyQBEjAQIACzACDX0kgbAJOhJfMK9wlwP4LsyIZgiMAITZ06FBLGzEMTbJHLLj6YEokEDUCFGBRqzHmN5YEpPDCMjvSmF5uYeOFoUbZ42UEgDXh5Hk0vjeSCe47hVhwrJkSCcSJAAVYnGqTZYkcAbfCSxYUHvKlAHv55Zflbm5DIEAhFgJ0JkkCESZAARbhymPWo0sgX+ElSz5nzhxdgGHNSIbwCVCIhV8HzAEJRIEABVgUaol5jA0Br4SXBPLOO+/oAuyjjz6Su7lVgACFmAKVwCzYIoDlvb755hvx66+/2jqfJ3lDgALMG46MhQSyEvBaeMnEPvvsM12AzZ8/X+7mViECFGIKVQazkpHA4MGDtedIr169Mh7nTn8IUID5w5WxkoBGAMJrwoQJIptx/bJly1zTwqLUWLewX79+AmkxqEuAQkzduklyzlauXKm/xMGedMGCBUnGEWjZKcACxc3EkkLAb+GVFI5xLCeFWBxrNZplOnz4sGjUqJFJgNWuXVvs3bs3mgWKWK4pwCJWYcyu2gRyCa+OHTuKfHq81C49c+eEwK5du+hHzAkwnuspAdh7devWzSS+5Izqtm3bimPHjnmaHiNLJ0ABls6Ee0jAMQEKL8fIeMG/CVCI8VYIg8Ctt96aUXxJEda1a1ca5ftcMRRgPgNm9PEmQOEV7/oNsnQUYkHSTm5av/zyi7j99tt18VW1alX9O8RXnTp19N/du3cXR44cSS4sn0tOAeYzYEYfTwIUXvGsVxVKlUuIDRw4kEscqVBREczDihUrtCXLZC9XpUqVNJMI+RtbOHeGHZjc17BhQ7Fo0aIIllb9LFOAqV9HzKFCBCi8FKqMmGeFQizmFRxg8b777jtx/fXXiyJFiujCqlq1amL16tXixx9/1PdBdK1du1asW7dOnH766fp+LIfWpUsXsWbNmgBzHf+kKMDiX8csoQcEKLw8gMgoXBGAEMM6n6VLl9b/EGXvREFBgWCPmCusibkIzy6ILXnPQEz17NlTYDYuQiYBhv2YCdm3b1+TaCtbtqzYt29fYtj5XVAKML8JM/5IE6DwinT1xSrzFGKxqs7ACoPZjJdddpkmpGDTBb9fxmAlwOQ5X3/9tYCD1mLFiomWLVuKQ4cOyUPc5kmAAixPgLw8ngSk8KpVq5b+5mh8g6Q7iXjWexRKRSEWhVpSK49w2Lxjx46MmcolwORFuO9okC9peLOlAPOGI2OJCQEIr4kTJwoKr5hUaIyLQSEW48oNsGh2BViAWUpMUhRgialqFjQbAQqvbHR4TGUCFGIq1476eaMAC6+OKMDCY8+UFSBgR3gtXbpUgZwyCySQnQCFWHY+PJqZAAVYZi5B7KUAC4Iy01COAAxTsw01dujQQVB4KVdtzJANAnaE2Pbt223ExFOSQIACLLxapgALjz1TDoEAhVcI0JlkKAQoxELBHrlEKcDCqzIKsPDYM+UACVB4BQibSSlFQAqxUqVKpc3olX7E2COmVJUFmhkKsEBxmxKjADPh4I+4EaDwiluNsjxuCVCIuSUX7+sowMKrXwqw8NgzZR8JUHj5CJdRR5oAhVikq8/zzFOAeY7UdoQUYLZR8cQoEKDwikItMY8qEMglxO6++27BoUkVasrfPFCA+cs3W+wUYNno8FhkCFB4RaaqmFHFCFCIKVYhAWeHAixg4IbkKMAMMPg1egQgvF588UVLz/V0JxG9OmWOwyFAIRYO97BTpQALrwYowMJjz5TzIEDhlQc8XkoCWQhAiA0ZMkRYzZrk0GQWeBE8RAEWXqVRgIXHnim7IEDh5QIaLyEBFwQoxFxAi+AlFGDhVRoFWHjsmbIDArmEV/v27cWSJUscxMhTSYAE7BCgELNDKbrnUICFV3cUYOGxZ8o2CEjhVbt27TQnkoULFxYUXjYg8hQS8IAAhZgHEBWMggIsvEqhAAuPPVPOQoDCKwscHiKBEAns3LnT0kasZMmSgjZiIVaOi6QpwFxA8+gSCjCPQDIabwhQeHnDkbGQgN8EKMT8JhxM/BRgwXDOlAoFWCYq3Bc4gWzCq1ChQhxqDLxGmCAJ2CNAIWaPk6pnUYCFVzMUYOGxZ8pCCAov3gYkEA8CFGLRrEcKsPDqjQIsPPaJTpnCK9HVz8LHmACFWLQqlwIsvPqiAAuPfSJThvCaNGmSyDSrkUONibwlWOiYEqAQi0bFUoCFV08UYOGxT1TKFF6Jqm4WlgR0ArmE2KBBg7jot04r+C8UYMEzlylSgEkS3PpCgMLLF6yMlAQiR4BCTM0qowALr14owMJjH+uUcwmvdu3a0XN9rO8AFo4EMhOgEMvMJay9FGBhkReCAiw89rFM2Y7w+uKLL2JZdhaKBEjAPgEIscGDB2dc9BsOXTk0aZ9lPmdSgOVDL79rKcDy48er/02Awou3AgmQgBsCdoTYjh073ETNa2wQoACzAcmnUyjAfAKblGgpvJJS0ywnCfhLgELMX75WsVOAWZHxfz8FmP+MY5kChVcsq5WFIoHQCVCIBVsFFGDB8jamRgFmpMHvOQlI4VWnTh0Bv12pHxjX08YrJ0aeQAIkkIMAhVgOQB4dpgDzCKSLaCjAXEBL4iUUXkmsdZaZBMInQCHmbx1QgPnLN1vsFGDZ6PCYtlYjPNezx4s3AwmQQJgEKMT8oU8B5g9XO7FSgNmhlMBz0OM1efJkCq8E1j2LTAIqE6AQ87Z2KMC85ekkNgowJ7QScC6FVwIqmUUkgRgQoBDzphIpwLzh6CYWCjA31GJ4TS7hdcUVV9C4Pob1ziKRQNQJUIjlV4MUYPnxy+dqCrB86MXgWgqvGFQii0ACJCByCbF77rlH0KFr+o1CAZbOJKg9FGBBkVYsHTvC6/PPP1cs18wOCZAACWQnQCGWnU/qUQqwVCLB/aYAC461EilReClRDcwECZCAzwQoxOwBpgCzx8mPsyjA/KCqYJwUXgpWCrNEAiTgOwEIMQw/YoHvVMfR2Dd06FCxe/du3/OhagIUYOHVDAVYeOwDSZnCKxDMTIQESEBxArD/shJipUuXTqwQowAL78alAAuPva8pHz9+PKsfL8xqpI2Xr1XAyEmABBQkQCFmrhQKMDOPIH9RgAVJO4C0KLwCgMwkSIAEIk9ACrGCgoK0ockk9YhRgIV3K1OAhcfe05QhvKZMmWLpuZ49Xp7iZmQkQAIxIbBt2zYxYMAAkVQhdvDgQa38YIDPrl27YlKz6heDAkz9OsqaQym86tatm/YWB4NTCq+s+HiQBEiABDQCQQgxGPu/9tpr4siRI6ROAoICLKI3AYVXRCuO2SYBElCagJ9CrGnTptqL8vDhw5VmwMwFQ4ACLBjOnqWSS3hdfvnlNK73jDYjIgESSCqBXEJs2LBhjt1XnHrqqZoA69mzZ1KxstwGAhRgBhgqf6XwUrl2mDcSIIG4EvBSiElfZLfccktccbFcDghQgDmAFcapdoTX4sWLw8ga0yQBEiCBxBDYunVrVmN9Oz1i5cuX13rArr766sRwY0GtCVCAWbMJ9QiFV6j4mTgJkAAJZCQAIda/f3/LWZPZhFiTJk00Ada4ceOMcXNnsghQgClW3xReilUIs0MCJEACGQi4EWI33XSTJsBOPPFEsW/fvgyxcleSCFCAKVLbFF6KVASzQQIkQAIOCEghVqJECU1cGdebhENXY48YfDXK49OnT3eQCk+NIwEKsJBrlcIr5Apg8iRAAiTgAQE7Quzbb78VRYsW1URYu3btPEiVUUSZAAVYyLX3xhtv6G9E8s0IW7iToHF9yJXD5EmABEjAIYFcQqx+/fraM/+EE04Qa9ascRg7T48TAQqwkGvzl19+EbJBUniFXBlMngRIgAQ8IpBNiMmX7S5duniUGqOJIgEKMAVq7dVXX2WPlwL1wCyQAAmQgNcEcgmxd955x+skGV9ECCgtwE4++WRtqm/x4sW1cfMiRYoIdNtiBgn2waldmTJlRNmyZUW5cuVExYoVRZUqVUT16tVFrVq1RI0aNbTfFSpUEKeccoqAQSQMJfFZtmxZRKqI2SSBcAn8+OOPWrtDm0N7Q1urVKmS1r7QxipXrqztO+mkk7S2hXb6zTffpGX6t99+0xaLR/vDufCJVK1aNVGzZk09DrRR2Mj06dMn7XruIIEoEsAox/z588Vtt92m/VfJ3i+5PeOMMwTaBkPyCCgrwI4dO5bRNkretPlsIeBWrlyZvNpmiUnABYH9+/eLZs2aaS8zhQsXtmyXeDnCovDnnnuu2LRpU1pKv/76q2jRooWoV6+egEjL1IYRP4bkR40alXY9d5BAlAj8/PPP4tlnn9U6AzLd63Lfgw8+SAEWpYr1MK/KCrC9e/fqD2j0aGHGCN6KR48eLV5++WUxc+ZMMWfOHPHuu++K999/P+Nn1qxZehzyZscWjYKBBEjAOYHdu3eLQYMGpbUr9IrhD8duOHLkiMCbv7Fd3nDDDWLHjh12o+B5JKAsgddee03rJTbe33jxhwPWXr16ibFjx2qTrDA8yZBcAsoKsM2bN+sP5+XLl7uqITzQjQ0A33v06OEqLl5EAiTwHwJdu3Y1tS30XGXq9frPFeZvn376qakXrHv37uwFMCPirwgSwFAivOQb/3cwTP/YY4+JXbt2+V6iQ4cOmcwFTjvtNK3nGiY5GOqvWrWqJgyxH6Y7GPIvKCjQhv3Rg43fyO/vf/978cc//lHcfvvt4pVXXhEbN270Pe9JTEBZAQYbLdzEsOtyMz7+3HPPmRoB4sIyELhBGUiABPIjALGV6niyX79+tiLdvn279kcg/6RatWoljh49autankQCKhNYsGCB/r8DG+ZnnnlGwJwmqACRd95552kCSq47KdtZvlsIs//+7/8WMElg8IaAsgIMw4u4YW688UbHJf344481Q33jDYchEqp4xyh5AQlYEhg8eLD+Z4O2VqxYsZxtDEOPLVu21K9r0KCB2LNnj2UaPEACUSKAFxNM+MLwOpyuhhnQcfHQQw/pbQ1tFL1c6NF6++23xbx587Rh0FWrVgmMMv39738XM2bMELBJu+KKK9L+Q+X/6amnnqqdG2bZ4pK2sgJswoQJ2o2DrZOwZcsWrddM3izYwuAXNxsDCZCAdwQwOxLDGMa21rt376wJ4IVKno/Zk3wpyoqLByNI4PDhwwIrnKgQYCst2xu2mCRjN/zrX/8S3bp1E5km3qD3m+4z7JK0Pk9ZAfbJJ5+I+++/X6xfv9469ylHMIyBWVbGGw7fR44cmXImf5IACXhB4KmnnjK1NysXFEgLMxtl24StiVvbTi/yzThIIAkEMAQq2xy2nTt3dlxstNPzzz/fFA/iwhArJ804xmm6QFkBZsqlzR/ws2K82fAdCp6BBOJC4KeffhJz584V8C2kQoBNJXx5Gdvdtddem5a19957Txv+wHmYDYbfDCQQZwIYAsQMfbTZsMJdd91lapsPPPCAq6xg9nMmmzI3JkKuMhDTi2IjwCZOnGi60fCgb9iwoTh48GBMq47FSiKBNm3aaPf5n//8Z2WKP3nyZFPbw5DFl19+qecPTlnxtixFGtoqAwnEncCkSZO0e/6ss85yNZHMCz6tW7fW2x3a31//+lfX0cK1hmzDcoseb5giMLgjEAsB9vnnn2ue8eVNgS0MIdetW+eOCq8iAUUJwGgd9zfsp1SZOYjeOPzJGNtf+/btNYIwsDf6+3L7Bq5odTBbJGBJAC9Jsk2g1zqMAHcTMg/Yrl69Oq9sdOjQwRQf4sRSegzuCERegO3cudM0pR03BGZ6oOuXgQTiRuDRRx/VH4Aq2TbOnj1bz5d84C9cuFBcdtll+n745WMggaQQQM+vNGBv3rx54GYDsM+SbRFbLPGVr0sMGN4b48R3OJZlcEcg0gIMM00uuuiitBvi8ccfd0eDV5GA4gR++OEHrXcXDz44UFTJkP3CCy9Ma4vyYX3JJZfk/fBXvGqYPRJIIwAbZNkG4BIiyICZ/zJt69vLhwAADBpJREFUbNFLnW9Yu3atKU7Ee+mll+YbbWKvj7QAGzBgQNrNcPXVVye2MlnwZBCYNm2aft9j8Xk8FFUIn332mZ4v44P/7LPPFlhajIEEkkZg27ZtuvE6esOCtH/EknvGdujFhDR0emASjTFeeM1ncEcgsgLM+CckbwbYxzhZj84uMiwizEACKhEYOHCg/hCEPRhWjlAhXHnllXq+ZLv0wxYTM8zcrJChAiPmIVkEMBRfvHhxrV3APCYo04Gbb77Z1Ba9mrhjtOlEG2/UqFGyKtTD0kZSgGGGFYZf5AMeW8yy8ronYMWKFQJT6uHhO8ypxB7WN6OKCQG8FGCdNtkGSpYsKaZOnRp66WB7KfMkt3gT9ypgVjOWGcO6dlgUnIEEokAAbleM/1lYk/jAgQO+Zj3VJ+asWbM8Sa9p06amNo4ebgZ3BCInwDDlFQuLyoc7tujaxdJFXgTM6EJcbdu2NaUB2xu/wocffijeeOMNfsjA0T0wffr0NE/VeGEIYtHfTG0Ba8TB9YuxbeI7li7Jd7mhzZs3i+HDh+vDOYi3b9++mbLhyT70pLNN8pnk5T2ABbmN7lhOP/10AYfjfgT0DsPZsbEtwrO9FwFrTRrj7dKlixfRJjKOSAkwvPXD4M9Y+fiOmWH5BqzbNXToUG16f2r8+O2nAGvcuHFamTLlgfsKkVOh3AyqV68euAsWvLjA9YTVPXrPPfc4bqJYN/L111/XZlJi6CY1bj8F2HfffZeWXmr6/J37XiSj3Iz86LnesGGD6f5FD7lXpjT16tUzxY01YRncEYiUABs2bJip4tG4O3Xq5IktiFxqAd3Ed955pzacaVznjgIs94OED1t1GL388svunggurzLapMH1BHqrjPcDhvGdLCuGbHzwwQd6HDD0/ctf/qLZz8h4KcDUud9knXDrvE6w8LXXAaM4xrrAsKFXwdiLhzTGjx/vVdSJiycyAgyrtBtvKHw/88wzxb59+zypNNipwEjRKLQqVqyop2nc70mChkjwh4XlXPghA6f3AGZBYmFcY9vAG2m+/n4Mt2fOr3gAy/RhHoC2giG8ChUq6PtxPNMSRdkixwwy+A7DdHoZxo0bp8fppwBDD4LTuuD5bL+57oFU8QLXEE5fTGRbyLZ94okn9HaCttezZ89sp9s+tmnTJlO8iNvYPm1HxBM1ApEQYP/85z/TxrPLlCkjvv76a1+rMSgB5mshGHlsCWCYzDgcULt2bd9sSqwgwn5RTkuHEDT6JUudBo+HNVxV5BOCEmD55JHXkkAqAbwQYd1EtAF80CP85JNPejYsmJren/70Jz0tpDdq1KjUU1z9njJliileDG1ygporlNpFygsw+A9KnfYKo3uvZnRkQ0cBlo0Oj4VJAAa16P2SD/R27drlbejutDx4ATIO07/00kumKGDDVaNGDT2PyCsMePMJFGD50OO1YRCA76yrrrpKbwdot4sXL/Y1K6l2xV6tDJMq7O644w5fyxH3yJUWYJjJ0bFjR/3GlX82Dz74YCD1QgEWCGYm4pAAZgLXrVtXbxe33npr4MucYJjRmIe77rorYylS35jRhjF7022gAHNLjteFRQBtQ/53wWwGM3r9DBB80u+YTHfLli15J4kXqvLly+tlQUcIllticE9AaQFmXPdO3kh403c7m+OWW24RMBa2GyjA7JLieUESgC2VbA8Y1gjaISkWATcuAYZeLSubM8yOrF+/vp5f5LtOnTquFxKnAAvyTmNa+RJ488039Xsfvuu2bt2ab5Q5r//qq6/0NNHeTjnllJzX2DnBaOuJeDHrmSE/AsoKMCz6CYUt/2iwhb2LW39CeANAfE6cxlGA5Xdz8WrvCXz88cd6mzjnnHNcC5l8cgaDXtku0UZy/anAl5I8X25Hjx7tKgsUYK6w8aIQCOClRPqshM1XUKtVoIdZtjNs//CHP+RdevSq4cXJGC+N7/PGKpQUYPBonzpbpFSpUmLNmjWuSzxmzBjt5nEyE4sCzDVuXugTAbm4L14mlixZ4lMq1tHCcFg+hGF8D0GYK6CHLtV7Nt7KMZTqNFCAOSXG88MiADtl2Vbc+MFzm++HH35YTxfpe2Gnldr7hdEkhvwJKCfAsDwDpubKG1du8RadT7jwwgu1OJ04baUAy4c4r/WDwGmnnabdx61atfIj+qxxYjjF2CvtZImhuXPnprVpJ+YAMmMUYJIEt6oTgJsU/H+hzeTqJfayLEaDf6Q/duzYvKJfuXKlaRklmB9YmRzklVACL1ZOgF1zzTVpD2p4qM8nYMaY/ONwIuQowPKhzmv9IDBx4kRxwQUX+O6CJTXvS5cuFZhyLl+InPQky7iks2MZR9GiRQVWoHASKMCc0OK5YRKA+6RmzZqJadOmBZqNVK8BCxYscJ0+TH6MQ48YUt29e7fr+HihmYBSAmzEiBH6A14+pDH0iGGORYsW2frMnz9fzJ49W7zyyisCa2/BkSPW3JLxORmHpwAz3yz8lTwCsP1AW0pdV+6ZZ55xDCOTXzDYdX7//fe246IAs42KJyaQwKFDh0Tqsl1unYijXcLOVP534hmwevXqBFL1r8jKCDD4RSlSpIhe2bLSvd7Cu7bdQAFmlxTPixsBLBKMGcew1bJqg1hkG0Oh2TzSw1akTZs22oM8VcTJeGFLBhsxLOqby6aMAixudxrL4yWBFStWmNprpUqVXEUPkwFcK9soZnB++umnruLiRdYElBFgcOIoK9uvLXyjOHFhQQFmfePwSLwJOGmPWPzbKmQTcJna+VtvvWUVlbafAiwrHh5MOAGsAWtsVxdffLFtIujtxsoWeKkyxnHTTTd5tuSf7cwk5ERlBBjWjtu4caOAu4idO3dqM6SwziOM8g8fPqwZ/cGnEAQUZlXhg98wBoRfIpyD7lecv3//fm15BIxV79ixQ4sTa1g5nXVFAZaQVsBiphGQ7RE9xhjCQJuCI0a0ObQ1rFCxa9cubfgwm2NJtGkMZaBNY8kStE+0VxkP2jjaKdLJFo/MIAWYJMEtCZgJoF21bdvWJJ5g+wy7rdatW4tevXqJRx55RMAUYPLkyZppAbwDYNZk586dxUknnWS6FhN+MPGGwT8Cyggw/4roPmYKMPfseCUJ+EGAAswPqowz6gSGDBkiKleubBJQxl4sJ9+xvNhtt92mdV5EnYvq+acAy1JDFSpU0G9ozvzIAoqHSCAgAi+88ILeJvv06RNQqkyGBNQmYPyvciK20EOG5YUaNmwounbtKmbOnKn1dKtd2vjkjgLMoi4x3ALvxfJmzscJrEUS3E0CJOCQgNHJJP4wGEiABIRmXoNhfAzzwxRHmufguzTJgcnA9u3bNZOADRs2aMP+sPtiCI8ABZiBPRYWhRHiq6++Kq688kpdfEGEwYfR1KlTBWaHwL8LAwmQgP8EDh48qM2MhGuZUaNGiXLlyuntEpNqYNMCOxX4OoLdGgMJkAAJRIUABZihpjANXvZ4ZdvC0JGBBEjAfwLwwp2tLRqPwZUNAwmQAAlEhQAFmKGm4PUXtl6Y4YU3b8wqwaxLzLREVy5mbGFGGI4xkAAJ+E8AQyRyBiWGUjADE/vQLmEmgF4vDLugXWJmNAMJkAAJRIUABVhUaor5JAESIAESIAESiA0BCrDYVCULQgIkQAIkQAIkEBUCFGBRqSnmkwRIgARIgARIIDYEKMBiU5UsCAmQAAmQAAmQQFQIUIBFpaaYTxIgARIgARIggdgQoACLTVWyICRAAiRAAiRAAlEhQAEWlZpiPkmABEiABEiABGJDgAIsNlXJgpAACZAACZAACUSFAAVYVGqK+SQBEiABEiABEogNAQqw2FQlC0ICJEACJEACJBAVAhRgUakp5pMESIAESIAESCA2BCjAYlOVLAgJkAAJkAAJkEBUCFCARaWmmE8SIAESIAESIIHYEKAAi01VsiAkQAIkQAIkQAJRIUABFpWaYj5JgARIgARIgARiQ4ACLDZVyYKQAAmQAAmQAAlEhQAFWFRqivkkARIgARIgARKIDQEKsNhUJQtCAiRAAiRAAiQQFQIUYFGpKeaTBEiABEiABEggNgQowGJTlSwICZAACZAACZBAVAhQgEWlpphPEiABEiABEiCB2BCgAItNVbIgJEACJEACJEACUSHwf0XfaHo8rztzAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.013931,
     "end_time": "2021-04-05T17:56:54.649609",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.635678",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The next graph is an example of J. Pearl (different notation), where the graph is considerably more complicated. We are interested in $D \\to Y$.\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "Here we try conditioning on $X_2$. This would block one backdoor path from $D$ to $Y$, but would open another path on which $X_2$ is a collider, so this shouldn't work. The application below gave a correct answer (after I put the spacings carefully).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.136442,
     "end_time": "2021-04-05T17:56:54.800269",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.663827",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "graph<- \"z1 -> x1\n",
    "z1 -> x2\n",
    "z2 -> x2\n",
    "z2 -> x3\n",
    "x2 -> d\n",
    "x2 -> y\n",
    "x3 -> y\n",
    "x1 -> d\n",
    "d -> m\n",
    "m -> y\n",
    "\"\n",
    "dosearch(\"p(y,d,x2)\", \"p(y|do(d))\", graph) #observed only (Y, D, X_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.013831,
     "end_time": "2021-04-05T17:56:54.828405",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.814574",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Intuitively, we should add more common causes. For example, adding $X_3$ and using $S = (X_2, X_3)$ should work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.239136,
     "end_time": "2021-04-05T17:56:55.081723",
     "exception": false,
     "start_time": "2021-04-05T17:56:54.842587",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "print(dosearch(\"p(y,d,x2,x3)\", \"p(y|do(d),x2, x3)\", graph)) #can ID conditional average effect?\n",
    "print(dosearch(\"p(y,d,x2,x3)\", \"p(y|do(d))\", graph)) #can ID unconditional effect?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.014749,
     "end_time": "2021-04-05T17:56:55.112369",
     "exception": false,
     "start_time": "2021-04-05T17:56:55.097620",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "This retrives correct formulas for counterfactual distributions of $Y(d)$ induced by $Do(D=d)$:\n",
    "\n",
    "The conditional distribution is identified by:\n",
    "$$\n",
    "p_{Y(d) \\mid X_2, X_3}(y) := p(y |x_2, x_3: do(d)) = p(y|x_2,x_3,d).\n",
    "$$\n",
    "\n",
    "The unconditional distribution is obtained by integration out $x_2$ and $x_3$:\n",
    "\n",
    "$$\n",
    "p_{Y(d) }(y) :=  p(y do(d)) = \\sum_{x2,x3}\\left(p(x_2,x_3)p(y|x_2,x_3,d)\\right).\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.014567,
     "end_time": "2021-04-05T17:56:55.141601",
     "exception": false,
     "start_time": "2021-04-05T17:56:55.127034",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Next we suppose that we observe only $(Y,D, M)$. Can we identify the effect $D \\to Y$?  Can we use back-door-criterion twice to get $D \\to M$ and $M \\to Y$ affect? Yes, that's called front-door criterion -- so we just need to remember only the back-door and the fact that we can use it iteratively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.281619,
     "end_time": "2021-04-05T17:56:55.438852",
     "exception": false,
     "start_time": "2021-04-05T17:56:55.157233",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "print(dosearch(\"p(y, d, m)\" , \"p(m|do(d))\", graph))\n",
    "print(dosearch(\"p(y, d, m)\" , \"p(y|do(m))\", graph))\n",
    "print(dosearch(\"p(y, d, m)\" , \"p(y|do(d))\", graph))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.016414,
     "end_time": "2021-04-05T17:56:55.472018",
     "exception": false,
     "start_time": "2021-04-05T17:56:55.455604",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "So we get identification results:\n",
    "First,\n",
    "$$\n",
    "p_{M(d)}(m)  := p(m|do(d)) = p(m|d).\n",
    "$$\n",
    "Second,\n",
    "$$\n",
    "p_{Y(m)}(y) := p(y|do(m)) = \\sum_{d}\\left(p(d)p(y|d,m)\\right), \n",
    "$$\n",
    "and the last by integrating the product of these two formulas:\n",
    "$$\n",
    "p_{Y(d)}(y) := p(y|do(d)) = \\sum_{m}\\left(p(m|d)\\sum_{d}\\left(p(d)p(y|d,m)\\right)\\right) \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.015846,
     "end_time": "2021-04-05T17:56:55.503891",
     "exception": false,
     "start_time": "2021-04-05T17:56:55.488045",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The package is very rich and allows identification analysis, when the data comes from multiple sources. Suppose we observe marginal distributions $(Y,D)$  and $(D,M)$ only. Can we identify the effect of $D \\to Y$. The answer is (guess) and the package correctly recovers it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.592959,
     "end_time": "2021-04-05T17:56:56.112631",
     "exception": false,
     "start_time": "2021-04-05T17:56:55.519672",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "graph<- \"z1 -> x1\n",
    "z1 -> x2\n",
    "z2 -> x2\n",
    "z2 -> x3\n",
    "x2 -> d\n",
    "x2 -> y\n",
    "x3 -> y\n",
    "x1 -> d\n",
    "d -> m\n",
    "m -> y\n",
    "\"\n",
    "print(dosearch(\"p(y,m)\\np(m,d)\", \"p(m|do(d))\" , graph))\n",
    "print(dosearch(\"p(y,m)\\np(m,d)\", \"p(y|do(m))\", graph))\n",
    "print(dosearch(\"p(y,m)\\np(m,d)\", \"p(y|do(d))\", graph))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.7"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 68.035072,
   "end_time": "2021-04-05T17:56:57.154022",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-04-05T17:55:49.118950",
   "version": "2.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
